Modeling of diffusion processes on physical and 
structural levels of materials 

A. L. Svistkovt A. V. Ilinykh 

Institute of Continuous Media Mechanics, 
Academic Korolev Str. 1, 614013 Perm, Russia 



Abstract 

A possibility to use an integral operator for establishing the link be- 
tween physical and structural levels of materials in modeling diffusion 
processes is considered. We show how to perform the transition from the 
stochastic description of motion of a system of points to the continuum 
description of diffusive profiles. 

1 Introduction 

Materials used industrially are complex multilevel systems. Analysis of these 
systems generally involves the consideration of their physical, structural and 
macroscopic levels (micro-, meso- and macrolevels) . As a rule, for each level we 
can formulate a special equation that allows modeling the processes taking place 
in the material: molecular motions and interactions on the physical level, super- 
molecular formations and processes in the vicinity of inclusions in composites 
on the structural level, and the state of the body as a whole on the macroscopic 
level. 

The considered processes can be studied separately on each level, which 
essentially simplifies their analysis. However, this possibility does not always 
exist. Such a modeling can be realized in the case when the points of higher 
level are the material regions that are representative (contain a great number of 
elements) and in homogeneous state. Under these conditions, the transition to 
state parameters of higher level can be performed by the averaging procedure. 
The studied volumes must be representative enough to exclude the fluctuation 
of parameters. The homogeneity of the state of these regions makes it possible 
to search the functional dependences on the higher level. 

There is another way of establishing the link between processes on different 
levels. It suggests the construction of the mathematical operator for creation 
of process images on the higher level instead of the use of the hypothesis of 
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representativeness and homogeneity of volumes. For composites, the transition 
from the structural level to the macroscopic one by the appropriate mathemat- 
ical operator has been discussed in |Q , Q . This idea can also be used to realize 
the transition from the physical level to the structural one in modeling diffusion 
phenomena. In the present paper, we consider the simplest case that leads to the 
Fick law. The analysis of this case allows us to compare the results of numerical 
modeling of diffusion of a great number of particles with the process images on 
the level of continuous medium, to follow the disappearance of fluctuations and 
to find the most suitable form of the mathematical operator. 

2 Modeling of diffusion of dissolved material con- 
stituents on physical level 

We consider diffusion of n material points (particles). Their motion in the 
Euclidean space reproduces the diffusion of dissolved material constituents. It 
is assumed that the change in the position of the material point in space can be 
defined by Ito stochastic equations 

dx„ = v n dt + b n dz n , ra = l, ..JV, (1) 

v n = v„(xi, ...,xjv), b n = 6„(xi, ...,xjv), 

where x„ is the radius- vector of the n th particle at time t, v„ is the vector 
function, which defines the determinate component of particle motions in the 
Euclidean space, z„ are independent Wiener's processes, and b n > is the 
scalar function characterizing the peculiarities of the chaotic motion of the n th 
particle. 

Let us clarify the obtained expressions. It is assumed, in the framework of 
this model, that particles move chaotically and, as a result, redistribute through- 
out the material, i. e., the mass transfer takes place. Constitutive equations are 
given as follows. The random motion can be represented as the most proba- 
ble tendency (velocity v„) and the chance deviations from it. The vector v„ 
is a mathematical expectation of the rate of change of the n th particle space 
position. The deviations are described in the analysis of independent Wiener 
processes. The scatter value for random walks is given by the function b n . Our 
purpose in this work is to determine, in the framework of used mathematical 
models, conditions for the transition from the stochastic description of motion of 
points to the contunuum description of mass transfer by the diffusion equation 
of a material constituent. 

To derive the dependences we are interested in, the notion of a phase space F 
should be used, which is a set of possible values of the radius- vector components 
of particles. The phase space clement dr is determined from 

dr — dx-^ dx^ dx^ ... dx dx^j dxjy , 

where x l n is the i th coordinate of the radius-vector of n th particle. To every 
system state there corresponds the phase space point r. From here on the 
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symbol ip will be used to denote the probability density of particles at time t in 
phase space r determined by the radius- vectors Xi, xjy 



ip = ip{t,xx, ...,Xjy). 
The temporal variation of this density is defined by the Fokker — Plank equation 

f + E V x „- (V vn) - if] V x 2 „ : ( 6„ 2 l) = 0, (2) 

where 
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I is a unit tensor, and are the basis vectors of a rectangular Cartesian coor- 
dinate system. It is clear that for phase space points at infinity the probability 
density ip and its derivatives are equal to zero. 

Numerical modeling of the motion of particles in space is limited to prescrib- 
ing the random walks of particles Ax n within the time interval At according to 
the Gauss law ip n with the probability density [Q which has the form 
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2 , 
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where symbols v % n are the components of vectors v n in a rectangular Cartesian 
coordinate system. 

3 Fluctuating and determinate mass densities of 
diffusing material constituent 

In order to introduce continuum notions, we should perform the averaging pro- 
cedure over the vicinity of the point in the Euclidean space. The simplest, yet 
not the best, way allowing the introduction of the mass density of a diffusing 
material constituent suggests the division of the number of particles in the vicin- 
ity of the studied point by the volume of this vicinity. This procedure gives the 
function dependent on the radius-vector of the Euclidean space, which pulsates 
as the radius-vector changes. The derivatives of this function with respect to 
coordinates can not be determined in this case. The region of averaging must 
be large enough to suppress these pulsations. The transition to the smoothed 
differential continuum notions only be possible following the corresponding hy- 
pothesis. In this section, we consider more formalized way of inserting the mass 
density of a diffusing constituent in our model. 

Let us use the function 77 to introduce the mesolevel notions. By this func- 
tion one can evaluate the contribution of the n th material point having the 
radius- vector x„ at time t to the state characteristics of the continuous medium 
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at the point in the Euclidean space having the radius-vector x. We require that 
the function 77 be continuous and have the continuous first derivative and the 
piece- wise continuous second derivative and become zero on moving apart from 
the point x to the distance exceeding the value a r 

77(x-x„) = ^(x-x„) • (x-x n ) > a r , (3) 

satisfy the normalization condition 
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II(x 1 i 1 + x 2 i 2 + x 3 i 3 ) dx 1 dx 2 dx 3 = 1, (4) 
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and the requirement of independence of the space orientation of basis vectors 
IlixHi + x 2 i 2 + x 3 i 3 ) = II ( Q ■ + x 2 i 2 + x 3 i 3 ) 



where x l arc the coordinates of a rectangular Cartesian coordinate system, and 
and Q is the arbitrary rotation tensor. The function 77 must be continuous and 
twice differentiated in order we can find the continuous twice differentiated pa- 
rameters of the medium. The equations must satisfy the invariance requirement, 
i. e., they must be independent of the orientation of basis vectors. Hereafter, 
the function 77 of the argument x-x„ will be denoted by the symbol 77„ 

77„ = 77 (x - x„). 

Introduce the notion of the fluctuating density p for the diffusing material 
constituent by the following equality: 

N 

p(t,x) = myelin, 

n=l 

where m is the mass of one diffusing particle. Physically, the calculation of the 
density p involves averaging the particles over the volume of the material taking 
into account their distance from the considered point in the Euclidean space. 
The value p is the continuous and twice differentiated function, as the function 
77 is continuous and twice differentiated. 

Introduce the notion of the density p of the diffusing constituent of the 
material by the expression 

N 

p = m I n n ip dr 

p n=l 

which is a mathematical expectation of the fluctuating density. Good agreement 
between p and p is only possible when the parameter a r ([}]) is sufficiently large 
and the distribution of particles in the spherical vicinity of the radius a r of the 
considered point in the Euclidean space is almost uniform. 
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4 Derivation of the integral diffusion law from 
stochastic equations of particle motion 



Consider the mathematical expression 
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dr = o, (5) 



Its validity follows from (g), as the integrand function (in square brackets) is 
equal to zero. Let us rearrange equation (^|) using three identities. The first can 
be derived based on the time independence of the function 77^ 



The second is a consequence of vanishing of the probability density of state ip 
at infinity of the phase space r . 
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The last is derived from the expression 



_ x fc=l 
^ v fe=i 

V ; — 1 



dr 



dr 



5 



and has the form 
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With these identities, equality (||) becomes 



d_ 

at 



N 



E^W 



fe=i 



>, n=l V-=l ' ^ ' 



The replacement of the operator V Xn by the operator —V for the function II 
Vx.f^fli) = V x „77„ = -Vi7„, 



fc=l 



gives 
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where 



V... = j> 



dx l 



V z ... = V V 



Taking the operations 9/c?t and V out of the integration sign, we obtain the 
expression in question 
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p n=l 

It can be used to derive the sufficient condition for transition to modeling the 
mass transfer based on the diffusion equation of the medium constituent. 



5 Sufficient conditions for transition to contin- 
uum notions 

The diffusion process can not be modeled using the notion of the chemical poten- 
tial in the framework of the proposed approach. Such modeling must be based 
on the analysis of the first and second thermodynamic laws, the accounting of 
the interaction forces between the particles in the formulation of constitutive 
equations, and so on. With stochastic equations (|l|), the transition is possible 
to the Fick diffusion law only. 

Define the notion of a diffusion coefficient D by the equality 

JV 



d = j J J2n n b n 2 ipdr. 



To describe the process on the meso- level, it is sufficient to satisfy two conditions: 
the vectors v„ must obey the equality 

N 

pw + pVD = ml ^2ll n v n ipdr, (7) 

n=l 

where v is the center-of-mass velocity of the body and the functional relation 
must exist between two values p and D 

D = D(p). (8) 

With these requirements satisfied, equation (Q) can be written in the form 

dp 



ni , V>v) = V- LDVp 

which is a well-known Fick law. The possibility to describe the process by the 
equation of continuous medium is checked, in particular problems, by conditions 
(0) and (||). It should be noted that we are dealing not with the averaged values, 
but with their mathematical expectations. In practice it means that these values 
should be calculated by averaging over space and random realizations. 
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6 Integral operator kernel of third order accu- 
racy of image reproducing 

For problem solution we should define a particular form of the integral operator 
kernel 77. We take the distance r between the radius- vectors x and x* as an 
argument of function 77 

77 = 77(r), r = y/(x-x,) ■ (x-x»). 

To reproduce the effective density of particles most precisely, it is reasonable to 
use the function 77 (r) based on the ideas of constructing the integral operator 
kernel described by the RKPM method [|, §. 

The kernel 77 is supposed to meet the following requirements. Let some 
function /3(x) be represented as the coordinate dependence of the third order 
and the fast oscillating term 7(x) 

3 . 

(3(x) = 6 + b- x + B:x®x+B:x(g)x(g)x + j(x), 

3 

where b, b, B, B are tensors of zero, first, second and third rank. The cubic 
dependence must remain unchanged for any space point x* under the transition 
to the image performed by the integral operator 

J II(r) ^6 + b- x + B:x<g)x+B:x®x®x^Gn/ 
l 

3 . 

= + b • x» + B : x* ® x*+ B : x* <g> x* <g> x*. (9) 

The symbol V stands for the entire volume of the Euclidean space. In the 
problems where the fast pulsations of the function 7(x) suppress each other 
under integral transformation, i. c, 

J II(r)~((x)dV « 0, 
v 

the smoothed integral image of the initial function /3(x) will be the slowly chang- 
ing cubic function 

J II (r) /3(x) dV w b + b ■ x* + B : x* ® x*+ B : x* ® x* ® x* . 

V 

In this case the integral operator will work as a "filter" that keeps slowly chang- 
ing regularities and cancels fast pulsations. 

Determine the conditions which the integral operator of the third order ac- 
curacy should satisfy. The function /3(x) may be developed as a scries in powers 
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of Ax about point x 

3 . 

/3(x) = 6 + b- x + B:x(g)x+B:x(g)X(g)x- 



where 



3 . 

- a + a • Ax + A : Ax ® Ax+ A : Ax ® Ax ® Ax + 7(x) 



Ax = x - x* = AxMi + Ax 2 i 2 + Aa; 3 i 3 , 



3 

a, a, A, A are tensors of zero, first, second and third rank. We note that 
equality (^|) will be true only when the reproducing conditions are fulfilled 

Jn(r)dV = l, Jn(r)AxdV = 0, (10) 

v v 

J n{r) Ax® AxdV = 0, J Z7(r) Ax ® Ax <g) Ax = 0. (11) 

It is easy to find that in the spherical system of coordinates 

Ax 1 = r cos(a) sin(0), 
Ax 2 = r sin(a) sin(0), 
Ax 3 = r cos(9) 



the reproducing conditions (|l0j), (11) are equivalent to the requirement that the 
following equalities are valid 

r 2 n(r) dr = — , y r 4 77(r) dr = 0. 
o o 
An example of the operator kernel of third order accuracy is the function 

in which the constants 70 71 are chosen such that equations (]Io|), ( O ) are 
satisfied and the symbol H(...) denotes the Heavyside function. 



7 Integral operator kernel of first order accu- 
racy of image reproducing 

The transition from physical to structural level may be performed using integral 
operators of first order accuracy. These operators must meet the requirement 
that only linear dependences retain in integral transformation 

J n(r) (^b + b-x^dV = 6 + b-x*. 

v 
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An example of integral operator kernels of first order accuracy is the function 

t3 \ ri 2 



and 



n = < 



„3 




COS 7T 




< r < au 



< r < a r 



(12) 



0, a r < r, 

in which constants 7 are determined from the normalization condition 



r 2 dr 



1 

47T 



8 Numerical experiments 

In this section, we consider the problem of penetration of diffusing particles into 
a semi-infinite material with constant concentration of particles on the material 
boundary. The constant b will be used as a function b n 

b n = b, n = l, ...,7V, 



and all vectors v and v„ are assumed to be zero 



v = 0, 



v„ = 0, 



71 = 1, ...,N. 



In this case, equations (0) and (||) are valid. The diffusion coefficient is found 
by the formula 

b 2 

The performed numerical experiments demonstrated that the best smoothing 
effect can be achieved with the aid of the function II (r) of first order accuracy of 
image reproducing ( |l2] ) . This function must be constant in the region occupying 
the half of the spherical volume of radius a r . In the present work, graphs (Fig. |l|, 
Fig. ^) have been plotted based on the relation (|H|). 

The obtained numerical results led us to the conclusion that the analytical 
solution of the Fick equation H is in good agreement with the fluctuating density 
p only in the case of a great number of particles in the averaging region. For 
example, about 5000 particles entered the radius sphere a r on the material 
boundary when we calculated the dependence shown in Fig. [jj.a. 

An essential divergence between the fluctuating density p and the analytical 
solution can be observed under random realization of the process, in the course 
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0.5 1 1.5 x/y/Dt 

b 



Figure 1: Comparison of the theoretical diffusive profile with the numerical 
results when the concentration of diffusing particles on the material boundary 
is (a) 5000 per unit volume and (b) 30 per unit volume 
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0.5 1 1.5 x/VDt 



Figure 2: Comparison of the theoretical diffusive profile with the numerical 
results when the concentration of diffusing particles on the material boundary 
is 30 per unit volume. The results were averaged over 30 random realizations 



of which 30 particles enter the sphere with radius a r on the material boundary 
(Fig. [jjb). However, one can also obtain the mesolevel images of the process. 
For this, it is sufficient to perform the averaging of the fluctuating density over 
thirty random realizations. The obtained density (p) agrees fairly well with the 
analytical solution of the Fick equation (Fig. ||) and is close to its most probable 
value (mathematical expectation, i. e. the density p). 



9 Conclusion 

The state parameters for the medium studied (mass density, diffusion coeffi- 
cients) can be calculated with the aid of operators acting on a multitude of 
microscopic parameters a%, apj and having the form 

N 

L\(ax,...,a N ) = n n a n . 

n=l 

and 

N 

L 2 (a 1 ,...,a N ) = / ^ n n a n ipdr, 

p n=l 

where a n is the characteristic of the n th paticle. The fist operator represents the 
averaging procedure in the vicinity of the considered point in the Euclidean space 
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and leads to the fluctuating twice differentiated variables. The second operator 
gives the determinate continuous twice-differentiated parameters, which can be 
used as the state characteristics of a continuous medium. Determination of 
diffusive profiles analogous to those on the mcsolevel based on the fluctuating 
mass density is possible provided that thousands of particles enter the region of 
averaging. However, there is another way of constructing the mesolevel images 
of the process. It consists in averaging over the fluctuating densities calculated 
in numerical experiments. The results achieved in this case agree fairly well 
with the calculation data obtained on the level of continuous medium with 
much smaller number of particles entering the averaging region. 

The work has been done under the FPP "Integration" (project 00-01, grant 
contest for students and young researchers of Perm region) . 
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